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Abstract 


A non-local boundary condition is formulated for acoustic waves in ducts without flow. The 
ducts are two-dimensional with constant area, but with variable impedance wall lining. Ex- 
tension of the formulation to three-dimensional and variable area ducts is straightforward in 
principle, but requires significantly more computation. The boundary condition simulates a 
non- reflecting wave field in an infi. it e duct. It is implemented by a constant matrix operator 
which is applied at the boundary of the computational domain. An efficient computational 
solution scheme is developed which allows calculations for high frequencies and long duct 
lengths. This computational solution utilizes the boundary condition to limit the computa- 
tional space while preserving the radiation boundary condition. The boundary condition is 
tested for several sources. It is demonstrated that the boundary condition can be applied 
close to the sound sources, rendering the computational domain small. Computational solu- 
tions with the new non-local boundary condition are shown to be consistent with the known 
solutions for nonreflecting wavefields in an infinite uniform duct. 
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1 Introduction 

The emerging field of computational aeroacoustics addresses the resolution of short wave- 
length and small amplitude motions of compressible fluids [1]. These problems are typically 
posed on unbounded domains, which present special problems for numerical techniques which 
use grids as there are a finite number of grid points to fit in an infinite domain. Typically 
the grid points are spaced far apart in the farfield, meaning significant loss of the accuracy 
necessary for computational aeroacoustics. 

To limit the number of gridpoints, the domain is separated into a bounded part, which is 
gridded, and an unbounded part, which is treated analytically. The approach of this paper 
is to assume chat the unbounded, or exterior, domain is governed by a wave equation with 
a known general solution in terms of an eigenfunction expansion. To couple the bounded, 
or interior, domain and the exterior domain, it is not necessary to completely specify the 
exterior solution at the interface. Rather, at the interface, a relationship between the modal 
components of the exterior solution variables is expressed via a matrix impedance operator. 
The end result is that a computational solution is performed in the interior domain with the 
exterior solution replaced completely by a special boundary condition at the interioi /exterior 
interface. In current parlance, the boundary condition is referred to as "nonlocal”, in con- 
trast to the classical "local’ Neuman or Dirichlet boundary conditions, or the well known 
approximations of pseudodifferential operators [2]. 

It is assumed here that efficient and accurate solutions to problems in computational 
aeroacoustics depend as much on the computational boundary conditions as they depend 
on the interior algorithms. If accurate boundary conditions are specified at the surface of a 
small interior, the computational domain then is usually small and the physical phenomena 
may be accurately captured with computational economy. This assumption is examined 
using several sources in an infinite duct. These problems have classical analytic solutions in 
terms of an infinite series of duct modes. 

I he boundary operator, the nodal impedance operator is shown to be a similarity trans- 
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formation of a diagonal matrix operator whose elements are the modal impedances of indi- 
vidual waves in the classical solution to the exterior operator. The transformation matrix 
is constructed by evaluating each mode in the exterior solution at each node in the interior 
solution. 

Section 2 of this paper defines basic equations used in the computation and formally states 
the non-local boundary condition. The derivation of the similarity transformation for the 
boundary operator is given in section 3. Results and discussion of several source solutions 
and frequencies are given in section 4. Appendix A describes the numerical method, a 
finite element procedure, and Appendix B gives the details of the infinite-duct eigenvalue 
computations. These eigenvalue computations are a crucial step, since the computation of 
the boundary operator requires the same number of eigenvalues as boundary node points 
This may be a large number for the high-frequency cases, so that provision must be made 
for evaluating an essentially unlimited number of eigenvalues. Conclusions, relative to the 
basic assumptions of the paper, are given in section 5. 
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Outgoing Wave Section 



Figure 1: Infinitely long three dimensional duct and coordinate system. 

2 Computational Solution 

Consider a three-dimensional rectangular duct filled with a homentropic fluid within the 
region ^ < y < j and = y~ < 2 < y as shown in figure 1. Acoustic waves in a duct are 
governed by the linearized equations of mass, momentum and energy [3]. The background 
pressure (p) and background density (p) 1 are assumed constant and the background flow 
velocity components are zero. Under these assumptions the sound speed c is related to the 
steady density and pressure by the equation 



The duct is assumed infinitely long in the axial direction with a known velocity source at 
the plane x - 0. The walls of the duct contain sound absorbing material whose material 
properties vary arbitrarily both around the circumference and along the axis of the duct for 
0 < x < L. Within the region L < x < 00 the material properties of the liner are assumed 
uniform, so that the wave field generated by the source will be that of an outgoing wave in 
'or equivalently “mean” or “ambient” pressure or density (see [3]) 
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straightforwaid, but requires significantly more computation. The two dimensional analysis 
to follow is obtained by assuming planar waves across the duct, so that the wave field in 
independent of z (terms with fall out of (3)). 

The elements of the two-dimensional counterpart of equation (3) are combined into a 
second order differential equation containing the acoustic pressure [3]: 


V 2 p + k 2 p = S 


( 8 ) 


where 


k = 


UJ 


(9) 


is the freespace wavenumber and V 2 is the Laplace operator in the (x,y) plane. A source S 
has been introduced in order to generate non trivial solutions. The wall boundary condition 
may be expressed in terms of the acoustic pressure only: 


Vp • n + ikpcfij) ~ 0 


( 10 ) 


in which V is the vector gradient, vector n the unit normal vector to the liner surface and 
is the vector dot product. 

There remains the duct termination condition at x = L. Because the portion of duct for 
x > L is considered infinite and uniform, the wave field in this sect ion of duct is expected to 
propagate out to infinity without reflecting. The boundary condition at x = L achieves this 
objective. The condition itself is a linear relationship between the acoustic pressure and the 
normal component of acoustic velocity at x — L of the form 

( 11 ) 
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Here {;>, } and {u ; } are M x 1 column matrices containing the acoustic pressure and the 
normal component of acoustic velocity, respectively, along the boundary nodes at x = L. A 
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method for determining the complex coefficients in the M x M matrix [Z tJ ] is given in the 
following section for the case of an infinite duct. The general form, however, is valid for any 
duct, and may include, for example, the effects of radiation conditions at the end of a finite 
duct [4]. 

Equations (8) and (10) are solved numerically using a finite element method. Details 
of the numerical implementation are given in appendix A. The results are used to evaluate 
the performance of the wall lining over the length L. The following expression is used to 
evaluate the acoustic intensity at a point (x,y) in the duct [7]: 


lU-s) = 


ip 


. dp" 
dx 


(13) 


where the superscript asterisk denotes the complex conjugate and 3? is the real part.. The 
attenuation of the lining in decibels is then obtained from 

10 l°g lo ~^ (14) 

where 

W(t) - f 2 h I(x,y)dy. ( 15 ) 

is the axial power. To obtain W (x), the integration in equation (15) was performed numer- 
ically from the computational solution. The accuracy of the boundary condition is assessed 
by comparing the numerically computed sound power to computed values of the sound power 
determined from truncated exact analytic series solutions of the reduced wave equation de- 
scribing outgoing waves in an infinite duct. 
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3 Infinite Duct Radiation Condition 


Modal solutions to the homogeneous form of equation (8) are 


Pm( X ,y) = L) 


( 16 ) 


where <p{y) is the mode function 


Qmiy) = Ae~ tk * y -f- Be +ik * v 


(17) 


The admittance boundary conditions define the eigen values k y and eigenvector [o b J T 
as solutions of a homogeneous matrix equation for the vector [a b J T . If the upper wall 
admittance is 3{H/ 2) = and the lower wall admittance is 3(—Hj 2) = ft, then from (3) 
along with boundary conditions (4)-(7) and (9) it follows that 


{k y - kpc9 u )e~' kyl1 ^ -{k v + *«,)e^] | a | | 0 | 


(18) 


. ~{k y + kpc3i)e +,k * H/ 2 {ky - kpcdi ) e ~' kyH/ 2 

The determinant of the above matrix must be zero, which is the eigenvalue equation 

t _ e -*kyH ( Vi / fcyj- kpc3A 

\ky + kpc0 u ) \ky + kpcdi ) 

A method of solving this transcendental equation for values k y is given in appendix B. The 
solution for the eigenvalues is a crucial step in the method. Fortunately, the solution method 
given in appendix B can easily generate thousands of eigenvalues (solutions of (19) ) so that 
eigenvalue generation is not a problem. 

Assuming a discrete set of eigenvalues kym, the corresponding modal functions 0 m (y) 
are orthogonal. The eigenvector [a b is normalized such that the inner product of two 
modes is a Kronecker delta function 

r+Hj 2 
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The pressure in a field of progressive waves in an infinite duct is 

p( x -'j)= ]T Pt<Prr{y)e~ ,k,m{Z ~ L) 
m=U 


(20) 


( 21 ) 


The corresponding axial velocity amplitude is 


u{x,y) = ]T u+0 m (y)e 




m=0 


P C m=0 V * 


( 22 ) 


If we define the modal impedance as 


— pe- 


rn 


and relate modal pressure and velocity amplitudes by 


Pm ~ 


(24) 


then the pressure and velocity fields cam be given in more compact forms as 

pi x >y) = t 2mui<Pm(y)e-' k ^ x - L) ( 25 ) 

m= 0 

u(x.y) = f. uUMr (26) 

mz.0 

Matrix expressions, assuming a finite number of modes M, for the pressure and axial velocity 
at x = L are 

% 

p(L,y) = l*m(y)J[2«H«£} (27) 

u(L,y) = L^m(y)J{«m} ^ 

Note that the modal impedance matrix [ ] is diagonal for the present case of an infinite 

duct. In the case of a finite duct, this matrix would contain the radiation impedances (4) ot 
the duct termination and would not be diagonal. This are the only differences between the 
infinite and finite duct cases. 

Now define the pressure and axial velocity at node points, or {/-grid lines i, as 

Pi = p(L,y t ) (29) 

u, = - .L-Vi) (3 °) 
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The pressure and axial velocity at each node point are then given in terms of the modal 
amplitude coefficients for the progressive waves in an infinite duct. 


{p,} = [<t> im ){Z m ]{u+ m ) 

(31) 

{^} = [ d’jm ] { fiin } 

(32) 

d’xm = ^441(3/1) 

(33) 


Eliminating the modal amplitudes gives an operator relation between the nodal values of 
pressure and axial velocity. 

{p,) = [Z„ Hi,} (34) 

where the nodal impedance matrix is defined as 

[Z, J ] = [<A, m ][Z m ][ 0 jm ]- 1 (35) 

Equation (35) defines the nodal impedance matrix used to construct the non-local boundary 
condition. This boundary condition insures continuity of the solutions at the interface be- 
tween the interior and the exterior domains. Furthermore, the boundary condition, although 
constant, insures continuity for all possible exterior solutions. In other words, it represents 
the general solution in the exterior domain. 


4 Numerical Results 


In this section we validate the accuracy of the non-local boundary condition. The precise 
construction is given in appendix A. Exact attenuations were computed from the exact series 
solution for outgoing waves developed ir appendix B. The exact series expression used 51 
modes in the series solution. Sample calculations are presented for both rigid and soft wall 
ducts and for several sound sources. Each source oscillates at 100. 500, 1,000 and 5,000 
Hertz. It is shown that the new condition can be applied very close to the source, rendering 
the computational domain small. All results were computed using a duct one meter in 
height with ambient values of [>,p and c. The computation uses a grid aspect ratio of 
unity (Ax - Ay) for each calculation with fifty-one evenly spaced points in the transverse 
direction. This grid was chosen because it allows approximately seven points per wavelength 
at the highest frequency of 5,000 Hertz. Although the calculations presented here will be of 
limited scope and represent only a small fraction of the total capability of the model, they 
do test the integrity of boundary condition formulation. The convergence characteristics and 
overall accuracy of the numerical model are also accessed. 

Each figure is presented in the form of the schematic of figure (2). The vertical dotted 
lines indicate the location of the right boundary where the nonlocal boundary condition is 
applied. The vertical scale is a measure of the attenuation from the axial x/ H = 0 to points 
x/H = AL/H,x/H — .2 L/H,...,x/H = 1.0 L/II consecutively. 

4.1 Results For A Rigid Wall Duct 

The first set of calculations is for a rigid wall duct (A = 0 + 0j). These choices for the 
wall admittance are not due to a limitation of the method, but allow comparison with 
exact analytical results available for planar and point sources. No attenuation is possible 
for a planar or point source propagating down a hard wall duct. Figure (3) compares the 
exact value of attenuation to that computed using the current model for a planar source. 
Numerical results were computed by applying the new boundary condition closer and closer 
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Figure 2: Axial locations where the nonreflecting boundary condition is applied, i.e. the six 
dotted vertical lines for x > 0 (at x/H = A, x/H = .2 etc.) are the points of application of 
the boundary condition. 
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Figure 3: Attenuation prediction for a planar source propagating in a hard wall duct. 
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Figure 4: Attenuation prediction for a point source propagating in a hard wall duct. 

to the sound source and computing attenuation over that length for each of the frequencies. 
The boundary condition was applied at the following values of duct aspect ratio ( L/H=l, 
L/H = .8, L/H=.6, L/H=.4, L/H = .'2 and L/H = .l ). Numerical predictions compare well to 
the exact value. There is ’ o attenuation in either case. 

Attenuation predictions for a point source in a hard wall duct located at y = 0 are com- 
pared to the exact series expression in figure (4). The point source attenuation comparisons 
are good for all frequencies and aspect ratios. Such good agreement for the point source 
is encouraging since a large number of planar modes are required for it full representation. 
Thus for rigid wall ducts, the boundary condition is accurate for a wide range of sources and 
may be applied close to the source. 
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Figure 5: Attenuation prediction for a lowest order mode source propagating in a lined duct 


4.2 Results With Wall Lining. 

A major goal of computational methods is to guide the design of liners which will reduce 
broadband noise in the current generation of aircraft. Any useful computational boundary 
condition, must remain valid in the presence of a wall lining. Attenuation comparisons are 
shown in figure (5) where the sound source is the lowest order mode. The admittance of each 
wall was chosen as /? = .5 — .5 i. Good comparison between the finite element and modal 
theory has been obtained in the presence of the wall lining. Note also that the presence of 
the liner causes attenuation of the sound as it propagates down the duct. 

Similar comparisons were obtained for a point source located at y — 0 with the same liner 

i 

\ values. Results for the point source are plotted in figure (6). Again, attenuation predictions 

using the nou-local boundary condition are generally in good agreement with results from the 
exact series solution. The exception is the 5000 Hz case where the computational solution 
does nol have monotone variation of the attenuation. Monotone variation is an absolute 
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Figure 6: Attenuation prediction for a point source propagating in a lined duct. 

requirement, since energy must be dissipated by the liner. It is posoible that the error 
observed here is due to the small number of points/ wavelength in this high frequency case. 
Further study will be needed to resolve this question. 














5 Concluding Remarks 

A non-local boundary condition has been developed for numerically simulating acoustic waves 
in ducts without flow. The non-local boundary condition represents the general solution to 
a linear wave equation in an unbounded exterior domain. The case of an infinite duct was 
used here for simplicity, but the boundary condition can represent the general solution for 
radiation from finite ducts as well. 

The boundary condition was applied to the case of two-dimensional, constant area ducts 
with constant impedance wall linings. Extension of the method to three-dimensional ducts 
is straightforward in principle, but requires significantly more computation. The boundary 
condition can be applied without change to variable-area ducts as long as the exterior domain 
represented by the boundary condition has constant area. Similarly, it can be applied to 
variable-impedance duct sections coupled to a constant-impedance duct termination. The 
boundary condition will represent finite ducts if the diagonal modal impedance matrix is 
replaced by a non-diagonal radiation impedance matrix. 

A finite element solution is developed which utilizes the boundary condition to limit 
the computational space while preserving the radiation boundary condition. The solution 
algorithm allows calculations for large numbers of computational points and high frequencies. 

The effectiveness of the new boundary condition has been evaluated by comparing pre- 
dicted sound attenuations with exact analytical results available from modal theory. The 
boundary condition was tested for both single mode and point sources at several frequencies 
using 3.4 or more points per axial wavelength. Excellent comparisons with exact analytical 
results were obtained with and without wall linings for most cases considered. The single 
exception was ihe highest frequency case in a lined duct which had only the minimum of 3.4 
points per wavelength. In this case, a small negative attenuation was observed in a part of 
the computational solution. 

It has been demonstrated that the boundary condition gives accurate results when the 
point of application is brought .Jose to the sound source. This is an important result in 
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computational aero-acoustics, since conserving grid points is a major concern. 

The boundary condition should be extended to the case of acoustic waves in the presence 
of steady flow. This extension will require the consideration of vortical and entropic modes 
which have been ignored here. The extension appears possible, however, since the only real 
requirement of the method is that the aeroacoustic field in the exterior domain is given by 
the general solution of a set of linear equations. 
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A The Numerical Method 


The numerical method chosen to solve enuation 8 coupled with the source, nonreflecting, 
and wall boundary conditions is a Galerkin finite-element method. Details on the method 
for structural vibrations problems are given in several texts [6, 5] and only sufficient detail 
is presented here to highlight important ditferences for the current duct acoustics problem. 
When applied to the current problem, the finite-element method may be interpreted as 
a physical visualization of the continuous pressure field as an assemblage of rectangular 
elements interconnected at nodes as illustrated in figure 7. It is assumed that there are N 
nodes in the axial and M nodes in the transverse direction of the duct as illustrated in the 
upper half of the figure. A typical rectangular element, [I,J] is shown in the lower half of 
the figure. Each element consists of four nodes labeled 1,2,3 and 4, respectively. A typical 
element is considered to have width a and height b as shown. The objective of the method is 
to obtain the unknown acoustic pressure at the nodes of each of the (M — 1 )(N — 1) elements. 

Galerkin’s finite element method is employed to minimize the error vector. This reduces 
the problem to a finite set of algebraic equations which are solved using matrix methods. 
Define the error function as: 

£(x,y) = V 2 p+* 2 p + S (36) 


Within each element p and S are represented as linear functions: 


P( J -y) = Ni(x,y)pi, S(x,y) = Y, W/(x,y)5/ 

i=i /= i 




iV<(x,y) = (1 - ~)~ 


(37) 

(38) 

(39) 


Ax 'Ay 

in which p/,5/, are the values of p and S respectively at node I. The variable admittance 
ll is represented in a similar manner along each boundary element. The correct solution to 
the sound field is obtained when the error E{x,y) is identically zero at each point of the 
domain. This is approximately achieved by requiring that the error function be orthogonal 
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to each basis function Ni(x,y) that are assumed to represent a complete set of functions. 
Contributions to the minimization of the error function from a typical element is 

L Jo E ( x 'y) N i( x ^y) d y dx - + {^/} ( 40 ) 

Where (.4^) is a 4x4 complex matrix (i.e. the stiffness matrix) {$'} is a 4x1 column vector 
containing the unknown acoustic pressure at the four nodes of the element, and {F, c } is a 
4x1 column vector containing source effects. The coefficients in the local stiffness matrix 
{A\j\ were computed in closed form and are not written explicitly in this work. The second 
derivative terms in equation (40) were integrated by parts in order that linear basis functions 
could be used and the effects of the wall boundary conditions included. 

Assembly of the global equations for the computational domain is a basic procedure in 
the finite element method. Appropriate shifting of rows and columns is all that is required 
to add the local element matrix [A' tJ } directly into the global matrix [A XJ \. A similar method 
is used to construct { } from {F'}. Assembling the elements for the entire domain results 
in a matrix equation of the form: 


[*«]{*,} = {*/} (4i) 

where [A] is an MN x MN complex matrix, and {<!>_,} is a A/A r x 1 column vector containing 
the values of the unknown acoustic pressures at the MN nodes of the duct. It is necessary 
to apply the source pressure, and termination condition to this system of equations before 
a unique solution can be obtained. Satisfying the noise source boundary condition consists 
simply of setting all nodal values of acoustic pressure at the source plane (j = 0) to the 
known value of source pressure. Provisions are made for insertion of this conditions into the 
assembled global matrix equation equation (41). Further details on this method of imposing 
source conditions are described elsewhere [8, 5, 6] and are not repeated here. 

1 he nonreflecting condition given by equation (11) must be imposed on the matrix equa- 
tion before the solution can be obtained. The axial velocity vector, {u} may be expressed 
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as a single equation on the acoustic pressure from the second of equation (3) 




(42) 


Substituting equation (42) into equation (11) gives a system of constraint equations relating 
the acoustic pressures along the nodes at x = L ( {p, } ) and those one grid point to the left 
of x = L ({Pjl}) 


[ZTMM = 


\M 

IU>p(xn — XjV-l) 


{Pjl} 


(43) 



1 

iup\XN — x^_i) 


[Z, } ] + [/] 


(44) 


where (/] is the identity matrix. Equations(43) are a set of constraint equations ensuring 
outgoing acoustic waves in the second section of duct. Standard finite element techniques 
are used to incorporate (43) into equation (41) using the Lagrange multiplier technique [5]. 

The global matrix [A,_,] generated by Galerkin’s Method following application of the 
source and constraint equations is a positive indefinite, complex matrix. Fortunately, owing 
to the discretization scheme used it will be block tridiagonal. The structure of matrix [A, ; ] 
is shown in figure 8. Note that it is a square block tridiagonal matrix whose order is 
MN . This global matrix contains a number of major blocks ([4/], [f?/] and [£>/)) which are 
themselves square and tridiagonal with the structure shown in figure 9. Note that each 
minor block [4/],[£?/] and [£>/] is a square tridiagonal matrix whose order is M. Much 
practical importance arises from this structure as it is convenient for minimizing storage 
and maximizing computational efficiency. Economy of storage is acheived by storing the 
rectangular array of coc !f icients within the bandwidth of [4 (J ] as shown in figure 10 

All computation,storage and boundary condition implementation is performed on this 
rectangular structure. Special matrix techniques exist for a solution of this structure. Gaus- 
sian elimination with partial pivoting and equivalent row infinity norm scaling is used to 
obtain the solution to this rectangular system. All computations were performed on Lang- 
ley’s Cray-2S computer. 
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| Figure 10: The stored MN x (2 M + 3) rectangular coefficient matrix with nonzero coefficients 
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B Eigenvalue Computation 

B.l High Order Eigenvalues 


When the walls are hard, the eigenvalues are multiples of x, but when the walls have finite 
admittances, a complex characteristic equation must be solved for the eigenvalues. This can 
easily be the most difficult part of the computation, so the method is developed here in 
considerable detail. We group the admittance at the wall with the frequency parameter kH 
by defining r 0 = pcfi 0 kH and r// = pc 3? . 

The eigenvalue computation procedures shown here were written independently of the 
rest of the paper and, purely by chance, used the time factor exp{-iu;<}. All complex 
variables with this convention are the complex conjugates of the variables in the body of this 
paper so that the eigenvalue equation 19 must be conjugated, giving the following eigenvalue 


equation. 


1 — e 2 ' kyH 


\ 1 + W V 1 + m7/ 

In the case of hard walls, where t 0 = t h = 0, solutions to equation 45 are 


1 ~ g jM / 1 ~ y7 

1 + v/j V 1 + 


k y H = mi, m = 0, 1,2, •• 


This elementary case suggests the definition 


(k y H) m - mx + 6 m , m = 0,1,2,--- 


B.1.1 Characteristic equation for higher modes, rn > 0 

The characteristic equation for the higher modes can be given in logarithmic form as 


F(6 m ) = 0 


F(f> m )=f>rn + ~logp + 


mx + 6 


T ° 1 . 1 i f, , T H 

~TT~ + o lo S 1 + — TT c 

T J ^ TTl *\ ~f" 


2 ,og [ 1 ~^m 


To 1 l 




It is clear from this form of the characteristic equation that the solution S m approaches zero 
in the limit where m oc. The points where rmr + fi m = ±r 0iH are singularities of F which 


ri 


cannot be solutions. These points should not cause difficulty except possibly in the case 
where m = 0. The explicit form of equation 49 is important. Computing each logarithim 
separately avoids crossing the branch cut, which could introduce an error of n in the value 
of the complex log. 


B.1.2 Newtons method solution for higher modes 


Newtons method numerical solutions of equation 48 utilize an iteration for S m given by the 
replacement operation 


■ 771 ' 771 


n 


F'(S m ) 


(50) 


where the derivative of the characteristic function is 

T 0 


TH 


F'(6 m ) = 1 - * 7 — T~ 1 ~ *7 T ( 51 ) 

(rrnr + 6 m ) 2 - r 0 J (m7r + 6 m y - r£ 

This solution method is very effective for higher modes m > 1. One can quickly generate 
a thousand or more eigenvalues on a modern workstation. The work in developing a useful 
basis for a modal computation is then limited to computing the lower eigenvalues. This 
computation is discussed below. 


B.2 Low Order Eigenvalues 

B.2.1 Alternate form of the characteristic equation 


The Newtons method is fast and dependable for the case where m > 1, but special care is 
needed for the cases where m — 0 or m = 1. The problem is worst when the admittances 
have large imaginary parts ! small real parts, so we attack these cases first. In these cases, 
there may be pure imaginary solutions for the eigenvalue K y H. Accordingly the following 
definitions are introduced to investigate these cases. 

K V H — -it] 
a = — t 0 t // 

. T 0 + T,/ 

-1 


b 


26 


2 


(52) 

(53) 

(54) 


I 


These changes of variable allow the characteristic eouation to be expressed as 

F(r)) = a + T) 2 + 2^/cotni; = 0 (55) 

To this point, there has been no real change, since al' variables are complex. Now consider 
solutions where the parameters a, 6 are real variables, that is, the real parts of the complex 
parameters defined above. We can find real solutions £ for these real parameters and add 
the effects of their imaginary parts later 

B.2.2 Solutions for real parameters 

Now take the real parts cq, 61 of the parameters a — oq + ia 2 , b = 61 -f i 6 2 , but, for brevity, 
use the same symbols a, b with the understanding that they are restricted to real values in 
this subsection. When the parameters are real, a better form of the characteristic equation 


utilizes the square of the eigenvalue. | 

f = , 2 =-(«„«)’ (56) 

f«) = « + ( + 2&| (57) 

C({) = cosh T) (53) 

sinh(n) 

5(0 = ^ (59) ' 

7 1 ' 


The value of the characteristic function at the origin and its behavior when £ is large are 
important for the determining if the characteristic function has zeros for positive 0 

F(0) = a + 26 (60) 

F{ 0 ~ 0 i-oo (61) 

The defined functions (7(0 and 5(0 have properties which make them useful in determining 
solutions to the characteristic equation. These propet ties will be listed here before consid- 
ering the numerical proceedure for solving for the eigenvalues. 
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ft 


The functions are defined by the series for the hyperbolic functions. 


0(0 

II 

» HM** 

(62) 

5(0 

S(2n+1)! 

(63) 


The series definitions permit the following immediate observations. 


5(0 

< 

0(0, 

0 < if < oo 

(64) 

5'(0 

< 

o'( 0 , 

0 < £ < oo 

(65) 

1 

< 

0(0. 

0 < f < OO 

(66) 

1 

2 

< 

o'( 0, 

0 < ( < oo 

(67) 

1 

< 

5(0, 

0 <( < oc 

(68) 

1 

a 

< 

5'(0, 

0 < £ < oo 

(69) 


The following formulas are useful for evaluating the derivatives of the functions. 

C(() = jS( 0 
~ m _ - £jj) 

Finally, the functions have the following asymptotic character. 


(70) 

(71) 


5(0 - C' /2 C( o, e-oo (72) 

no - \r l c( a e-oo (73) 

no ~ ^r ,/a c(o, *-»«> (74) 

The first derivative of the characteristic function is 


™ ■ -‘APPO 

(75) 

n o) = i + |* 

(76) 

no ~ i + b c x/t < ( - oo 

(77) 
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The derivative F'(£) approaches unity from above or below, depending on the sign of 6, 
wnen £ is large. The derivative F'( 0) may be positive or negative, depending on the value 
of 6. It is positive for b > —3/2, zero for b = —3/2, and negative for b < —3/2. The second 
derivative of the characteristic function is 


no 

F"{ 0) 

no 


b 


s 3 (0 



' C(Q(S 8 (Q- 2) + 5(Q 

2£ 2 


, (-*oo 


(78) 

(79) 

(80) 


The function within brackets is positive, so that the sign of the ond derivative depends 
only on the sign of b. If b is negative the curvature is positive, and if 6 is positive, the 
curvature is negative. The limit of the function within brackets is 4/45 for £ = 0, giving the 
value of the derivative at the origin. 

Small solutions |£| << 1 are possible for certain combinations of the paramaters u, 6. 
These combinations are identified by utilizing a two-term series for the characteristic func- 
tion. 


F( 0 = F(0) + F'(0)£ = 0 

which gives the following estimate for r/ 2 

‘-cm 


(81) 


(82) 


In order for this estimate to be consistent with the assemed smallness of £, the parameters 
must be related by 


- |F'(0)| e < F(0) < +|F'(0)|« 
e < F(0) < + 


•4 


*41 


(83) 

(84) 


In the special case where b — > -3/2, the limits above show that a — ► -26. These solutions 
may be positive or negative since they represent the effect of the parameters in moving the 
solutions from the real A V H axis to the imaginary h\H axis in the neighborhood of the 
origin. 
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There are a long list of cases to investigate in determining whether solutions for K y H 
exist near the origin or on the imaginary axis (£ is positive). The cases will be organized here 
by defining ranges of the characteristic function and its slope at the origin. These parameters 
( F(0) and F'( 0) ) depend on the parameters a, b only. 

1. If F'( 0) > -ft, the slope of the characteristic function is positive everywhere and, if 
F( 0) > +e, there is no positive solution, however 

(a) if|F(0)| < c, there is a small solution, £ = 0(e) , or 

(b) else if F(0) < — c, there is a positive solution £ > 0 or 

2. else if jF'(0)| < e,the slope at the origin is near zero and if F( 0) > +c 2 , there is no 

positive solution, however 

(a) if |F(0)| < e 2 , there is a small double solution £ O(e), or 

(b) else if F(0) < — e 2 , there is a positive solution, or 

3. else if F'(0) < — e, the slope at the origin is negative and the characteristic function 

has a minimum for some £ > 0. Find the location of the minimum £ mtn and evaluate 
the characteristic function to get its minimum value F min . If F min > c 2 , there is no 
solution, however 

(a) if |F m .„| < c 2 , there is a double solution near £ mm , or 

(b) else if F rniri < — ( 2 , the solutions depend on the value of F at the origin. 

i. If F(0) > +c 2 , there are two distinct solutions £ < £ min and £ > £ min , or 

ii. else if |F(0)| < c 2 , there is a solution near £ = 0 and another solution £ > £ min , 
or 

iii. else if F(0) <r - 1 2 . there is a single solution £ > £ m ,„. 
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B.2.3 Solutions for complex parameters 


Having found all possible solutions for real parameters, the next step is to find the effect 
of the imaginary parts of the parameters on these solutions. Let the parameters a , b be 
complex functions of a real parameter s which varies from zero to unity, 0 < s < 1 . The 
parameters are identical to their real parts when s = 0 and take their actual values when 
s=l. Begin with the solutions for real parts of the parameters in the case where s = U and 
trace these solutions to the point where s = 1. Note that these solutions can include effects 
of the real parts of the admittances in the parameter a. The tracing follows a differential 
equation for the square of the eigenvalue £(s) as a function of the parameter s. 


Ky(s)H = 


(85) 

a -- 

~t 0 th 

(86) 

b = 

■ T 0 + T H 
1 2 

(87) 

a(,s) = 

3?[a] 4- t§[a].s = a\ + ia 2 s 

(88) 

b(s) = 

+ i3[6Js = + ib 2 s 

(89) 


Now, £(s), 0 ( 5 ), and b(s) are complex functions of the real variable s. The characteristic 
function depends on s, and is complex, but is otherwise unchanged from the one defined for 
real variables. 



s) 

dF(i: 

f) 

e( 


t)'F(C; 

*) 

d? 


dF(i ; 

s) 

ds 


= <.(*) + { + 

(ClOSIO-l 

- 1+6(s) l-^— 


SH() 


C(0(S 2 (t) - 2) + S(0 


2V 


= a'(s) + 26'( .s ) 


C(Q 

S(() 


(90) 

(91) 

(92) 
(931 


The solutions where .s = 0 are initial conditions for a first-order differential equation which 
is used to find £(.s). 

(94) 


di = -i 4^rds 


m 
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This equation is valid as long as the initial condition is not a double eigenvalue. When the 
initial condition is a double eigenvalue £ m , the partial derivative in the denominator is 

dF d*F({;s) l 


nr 


d ( 2 


■U-U) 


(95) 


The differential equation for £ is then singular, but a differential equation for the square of 
the displacements of the eigenvalues from the double eigenvalue is not. 

(—) 

d(( - (mV = / a i F \ds (96) 

Va< J ) 

Given an initial step As, the two eigenvalues near the double eigenvalue are 





(97) 


The differential equation for the eigenvalue is then esed to trace each of these to the final 
value for s = 1. 


32 


REPORT DOCUMENTATION PAGE 

form Approved 
OMB No 0704-0188 

'^ocrt euroen *cr :rm :■ eerie" * rt’crrutior s ^st • 3 ;>e f '**^oc r 'sp. **c ri : 'cr -e. ^A n.; r»«r-aic"y veer^** rq ** aju *oufL« 

$oa 'n<*irt3ir>*ng the Odt<* ^eeaea. jna •:cn*Die'if' ■] irg re^e*."'.; :re . :iif!: ■: f '"-.it': 1 ' )*»na „ r-r^erts re^a'Si^q *h s Durden Mtirrate :r in. ;ther isoect ?f 

v ji.ea-on j» n for manor*. nciuainq suggestion* for reout m<g trvs ouraen ;■; Aasr.r von ^aac-jr'e's >*rvicev . r^cion te ’’or -* ; rmat cr* OoeMtion* tr-o 8 eoorT> U * 5 jetterson 
Cans Highway. »aite ’2C4. Arlington. .i ara to the O**" e ; f Management tra d^age: *> tperwcr* ^ecuct .c* ^'Cie-ct ;3.'C4 3 56) .Vasningt:r\ "C *3503 

1. AGENCY USE ONLY (Leave OUnk) 2. REPORT DATE 3. REPORT TYPE AND OATES COVERED 

March 1994 Technical Memorandum 

4. TITLE AND SUBTITLE 

A Non-local Computational Boundary Condition for Duct Acoustics 

5. FUNDING NUMBERS 
505-59-52 

6. AUTHORfS) ‘ 

William. E. ZorumsKi 
Willie a Watson 
Steve L. Hodge 

7. PERFORMING ORGANIZATION NAME(S) AND ADORESS(ES) 

NASA Langley Research Center 
Hampton, VA 23681-0001 

8. PERFORMING ORGANIZATION 
REPORT NUMBER 

9. SPONSORING /MONITORING AGENCY NAME(S) AND AODRESS(ES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 

10. SPONSORING MONITORING 
AGENCY REPORT NUMBER 

NASA TM- 109091 

11. SUPPLEMENTARY NOTES 

12a. DISTRIBUTION AVAHA3IL.TY STATEMENT 1 

Unclassified-Unlimited 
Subject Category 71 

1 

12b DISTRIBUTION CODE 

13. ABSTRACT (Maximum ^C0 words} 

A non-local boundary condition is formulated for acoustic waves in ducts without flow. The ducts are two-dimensional with 
constant area, but with variable impedance wall lining. Extension of the formulation to three-dimensional and variable area 
ducts is straightforward in principle, but requires significantly more computation. The boundary condition simulates a 
non-reflecting wavefield ir an infinite duct. It is implemented by c constant matrix operator which is applied at the boundary 
of the computational domain. An efficient computational solution scheme is developed which allows calculations for high 
frequencies and long duct lengths. This computational solution utilizes the boundary condition to limit the computational 
space while preserving the radiation boundary condition. The boundary condition is tested for several sources. It is 
demonstrated that the boundary condition can be applied close to the sound sources, rendering the computational domain 
small. Computational solutions with the new non-local boundary condition are shown to be consistent with the Known 
solutions for nonreflecting wave! lekJs in an infinite uniform duct. 

14. SUBJECT TERMS 

Computational Methods, Computational Boundary Conditions, Computational Aercacoustics 

15. NUMBER OF PAGES 

34 

16. PRICE COOC 

A03 

17. SECURITY CLASSIFICATION 

OF REPORT 

Unclassified 

18 SECURITY CLASSIFICATION 
OF THIS RAGE 

Unclassified 

19. SECURITY CLASSIFICATION 
OF ABSTRACT 

Unclassified 

20. LIMITATION OP ABSTRACT 


NSN 7540 0* 7*0 SS00 Standard Form 298 (R«v 2 89) 


bv -»NM Sid /If >* 

•m • • • 


*•* « » 



